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1.  INTRODUCTION 


This  report  describes  a  computer  program  developed  as  part 
of  the  contract  F33615-78-C-3412 ,  Trunk  Flutter  Analysis.  The 
computer  program  simulates  behavior  of  a  two-dimensional  trunk 
segment  in  presence  of  air  flows  existing  in  an  air  cushion 
landing  system  (ACLS) .  Through  such  simulations  a  greater 
understanding  of  trunk  flutter  mechanisms  can  be  achieved  and 
ways  to  eliminate  flutter  in  future  designs  can  be  developed. 

The  trunk- fluid  flow  system  addressed  by  the  program  is 
shown  in  Figure  1.*  As  shown  in  the  figure  the  trunk  is  assumed 
to  be  fed  by  a  fan  with  the  user  selected  characteristics.  The 
air  supplied  to  the  trunk  by  the  fan  flows  out  at  two  places. 

a.  Through  the  trim  valves  to  the  cushion 

b.  Through  orifices  at  the  side  and  the  bottom  of  the  trunk. 

The  cushion  air  flows  to  the  atmosphere  from  the  bottom  of 
the  trunk . 

The  trunk,  assumed  to  be  made  of  an  elastic  membrane  with 
a  finite  mass  and  flexural  rigidity  is  divided  into  a  number  of 
mass  nodes,  each  connected  to  the  other  by  springs  as  shown  in 
Figure  2. 

The  separation  point  is  assumed  to  be  always  at  a  particular 
slope  of  the  trunk.  However,  in  view  of  the  additional  work  that 
needs  to  be  done  in  identifying  the  location  of  the  separation 
point,  the  computer  program  provides  various  options  in  defining 
the  separation  point  location. 


*  See  Final  Report  "Contract  No."  F33615-78-C-3412  for  details  of  the  model. 


Figure  1.  ACLS  trunk 
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For  this  configuration  of  the  trunk-airflow  interaction  the 
computer  program  simulates  the  trunk  behavior  for  a  variety  of 
fluid  flow  and  trunk  parameters.  A  typical  output  of  the  simu¬ 
lation  consists  of  the  shape  of  the  trunk  defined  by  the  location 
of  the  nodes  as  a  function  of  time.  In  addition,  the  fluid  param¬ 
eters,  such  as  the  cushion  and  the  trunk  pressures,  and  various 

flows  are  also  available  |n  the  output. 

« 

The  computer  program  is  designed  to  be  flexible  and  versatile 
so  that  the  trunk  behavior  can  be  studied  for  variations  in  a 
number  of  parameters.  This  way  the  program  can  be  of  assistance 
in  developing  flutter  free  configurations.  A  list  of  major 
parameters  that  can  be  varied  is  shown  in  Table  1.  In  order  to 
ensure  that  the  user  can  employ  this  flexibility  easily,  the 
program  is  organized  in  such  a  manner  that: 

a.  The  major  parameters  are  varied  just  through  data 
input 

b.  Any  changes  in  the  program  made  necessary  due  to 
additional  experimental  information  on  flutter  char¬ 
acteristics  could  easily  be  performed. 

The  modular  structure  of  the  program  which  makes  such 
changes  easy  to  accomplish  is  described  in  Section  2  of  the 
report.  The  instructions  for  using  the  program  are  in  Section  3 
whereas  Section  4  illustrates  the  program  capability  through 
an  example.  Appendix  A  summarizes  the  equations  incorporated 
in  the  program.  These  equations  are  based  on  the  model  described 
in  the  final  report.* 


k r;<  e  Filial  Report  cl  Contract  No.  F33f;15-7S-C-3412  for  details. 


TABLE  1.  THE  CAPABILITIES  OF  THE  TRUNK  FLUTTER 
SIMULATION  PROGRAM 


The  program  simulates  the  trunk  behavior  for  variations  in  the  following  parameters 


ACLS  Parameters 

Trunk  Parameters 

attachment  points 
cross  section  length 

elasticity  variations  along  the  length 
flexural  stiffness  variations  along  the  length 
density  variations  along  the  length 
trim  valve  size 

trunk  orifice  size  and  location 
Fluid  Parameters 

fan  characteristics 
cushion  volume 
trunk  volume 
separation  point 
global  damping 

Operation  Parameters 

hard  surface  clearance 
Flutter  Suppression  Parameters 

strake  at  any  location 

minimum  gap  area 

external  spring  at  any  location 

Program  Parameters 

Simulation  Parameters 

time  step 
time  limit 
plotting  options 

Trunk  Model  Parameters 
number  of  nodes 
Options 

separation  point  options* 

a.  diffuser  model,  i.e.  separation  at  a  fixed  slope 

b.  fixed  gap  to  separation  point  height 

c.  trunk  orifice  flow  induced  separation  i.e., 
separation  occurs  at  the  last  orifice  row  if 

it  is  at  a  slope  less  than  the  diffuser  model  slope 

cushion  -  trunk  pressures  options 

a.  fixed  cushion  and  trunk  pressures** 

b.  variable  cushion  pressure  fixed  trunk  pressure** 

c.  variable  cushion  and  trunk  pressures  with 
fan  characteristics 

pressure  profile  on  trunk  bottom  options 

a.  nonvariable  pressure  profile** 

b.  variable  pressure  profile  without  trunk 
orifice  flow** 

c.  variable  pressure  profile  with  trunk 
orifice  flow 


*  For  flutter  studies  conducted  before  further  investigations  in  the  separation 
point  location  are  performed 

**  Used  for  initial  studies  on  the  behavior  of  a  particular  trunk  design 


Appendix  B  describes  a  program,  which  was  also  developed 
as  a  part  of  this  contract,  to  calculate  the  eigenvalues  of  a 
linearized  trunk  model.  This  program  can  enhance  the  under¬ 
standing  of  the  dynamics  of  the  trunk  motion  through  prediction 
of  the  natural  frequencies  and  the  mode  shape. 

Appendix  C  summarizes  Principal  Program  Nomenclature,  and 
Appendix  D  has  the  listings  of  the  computer  programs. 


2.  PROGRAM  ORGANIZATION 


Overall  structure  of  the  computer  program  developed  for 
simulating  the  dynamic  behavior  of  the  trunk  of  air  cushion 
landing  system  (ACLS)  is  described  in  this  section.  Details  of 
the  eigenvalue  computer  program  are,  however,  described  in 
Appendix  B. 

The  computer  program  has  a  modular  structure,  that  is, 
there  is  a  main  program  which  coordinates  operations  of  a  number 
of  subroutines,  each  of  which  perform  a  specific  function.  Such 
a  structure  makes  the  program  efficient  and  easy  to  modify.  As 
shown  in  Figure  3,  there  are  four  steps  in  the  program  execution: 

a.  Data  input  and  initialization 

b.  Initial  condition  acquisition  and  estimation 

c.  Dynamic  part  execution 

d.  Plotting  the  results. 

Details  of  each  of  these  steps  are  described  in  the  following 
subsection. 

2.1  Program  Execution  Steps 

The  manner  in  which  the  main  program  executes  each  of  the 
above  four  steps  is  shown  in  Figure  4.  The  main  program,  DYSYS, 
controls  the  dyn-.mic  simulation  of  the  trunk  model.  DYSYS 
coordinates  integration  of  the  differential  equations,  printing 
of  the  state  variables,  and  plotting  the  results.  It  initially 
calls  subroutine  EQSIM  to  initialize  values  of  the  derivatives 
of  the  state  variables.  DYSYS  prints  out  the  initial  conditions 
and  then  enters  the  integration  loop  which  calls  subroutine 
RKDIF .  DYSYS  calls  RKDIF  at  every  time  step  until  the  time 
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Figure  4.  Simulation  flowchart. 
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limit  of  the  simulation  is  reached.  RKDIF  is  the  integration 
subroutine  which  incorporates  a  fourth  order  Runge-Kutta  scheme. 
The  integration  scheme  requires  updating  the  differential  values 
of  the  variables  four  times  every  time  step,  therefore,  sub¬ 
routine  EQSIM  is  called  four  times  by  subroutine  RKDIF. 

RKDIF  is  the  numerical  integration  subroutine  which  calcu¬ 
lates  the  values  of  the  state  variables  at  time  t  +  dt,  given 
the  values  at  time  t,  using  a  fourth  order  Runge-Kutta  method. 
The  integration  scheme  is  summarized  below: 

a.  The  iteration  procedure  starts  with  the  values  of 
the  state  variables  y^,  y2  etc.,  at  time  t. 

Y.(t) ,  I  =  1 ,  n 

b.  The  slopes  Dy^U)  are  then  determined  from  y ^  ( t )  by 
calling  EQSIM 

Dy i  ( t )  =  dy ^ ( t) /dt 

c.  The  values  y ^  at  time  t  +  dt/2  are  then  determined, 

yil  =  yi  +  Dyi  ‘  dt/2 

d.  The  slopes  Dy^  (t  +  dt/2)  are  then  determined  by 

calling  EQSIM  and  using  the  values  of  y^  found  in  c. 
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e.  The  values  at  time  t  +  dt/2  are  then  determined 

yi2  =  yi  +  Dyn  ‘  dt/2 

f.  The  slopes  Dy^  (t  +  dt/2)  are  then  determined  from 

EQSIM  using  the  values  of  y ^  found  in  e.  above. 

g.  The  values  y^3  at  time  t  +  dt  are  then  determined, 

yi3  ■  yi  +  °yi2  •  at 

h.  The  slopes  Dy • 3  at  time  t  +  dt  are  then  determined 
from  EQSIM  using  the  values  of  y^  found  in  g.  above. 

i.  Finally,  the  values  of  the  state  variables  at  time 
t  +  dt  are  found  as  follows: 

yi(t  +  dt)  =  yi(t)  +  (Dyi  +  2Dyil 
+  2Dyi2  +  Dyi3)  dt/6 

During  each  integration  step  (that  is,  to  advance  from  t 
to  t  +  dt) ,  EQSIM  is  called  four  times  to  determine  the  slopes 
(b,  d,  f ,  and  h  above) . 


Subroutine  EQSIM,  which  calculates  the  various  flows, 
pressures,  motions,  and  forces,  is  the  only  model  specific 


subroutine  in  the  simulator  program  and  has  all  the  parameters, 
variable  initialization  input  and  system  equations  contained  in 
it.  This  subroutine  calls  subroutine  TRUNK  in  order  to  calculate 
the  initial  trunk  shape  for  the  given  trunk  length  and  attachment 
points . 

Once  the  simulation  is  completed  the  main  program  calls 
subroutines  PLOTTER,  PACKER,  PRNTPT  and  PSTORE  in  order  to  pro¬ 
duce  time  history  plots  of  any  of  the  system  state  variables  on 
a  printer  plot.  Table  2  summarizes  the  subroutines  used  in  the 
program. 


TABLE  2.  A  SUMMARY  OF  SUBROUTINES 


No. 

Subroutine 

Primary  Function 

Group 

1 

DYSYS* 

Main  program;  control  integration 
of  state  equations  and  I/O 

Main* 

2 

EQSIM 

Compute  state  derivatives  and  system 
pressure-f low-geometry 

Dynamic 

3 

RKDIF 

Coordinates  Runge-Kutta  integration 
algorithm 

Dynamic 

4 

PLOTTER 

Controls  data  plotting 

Plot 

5 

PACKER 

Read  plot  data  from  simulation 
output  file 

Plot 

6 

PRNTPT 

Plot  data  on  printer 

Plot 

7 

PSTORE 

Write  simulation  output  file 

Plot 

8 

TRUNK 

Compute  TRUNK  shape 

Geometry 

★ 

DYSYS  is  the  main 

calling  program. 
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3.  PROGRAM  USER  INSTRUCTIONS 


3.1  Program  Input  Data 

The  required  program  input  data  are  supplied  to  the  program 
in  three  ways: 

a.  By  data  cards  for  parameters  and  design  variables 
which  are  frequently  changed. 

b.  By  data  specifications  included  within  the  program, 
for  example,  physical  constants. 

c.  Through  subroutine  TRUNK  which  can  be  used  to  compute 
a  trunk  shape. 

3.1.1  Input  Data  Format 

The  input  data  set  format  specification  is  designed  to  allow 
flexibility  in  the  program  parameter  initialization  and  option 
specification.  Program  variables  are  input  in  a  specified  format 
and  sequence  using  FORTRAN  formatted  I/O.  Some  input  variables 
are  not  required  for  certain  operating  conditions.  These  special 
cases  are  noted  in  the  input  description  (marked  with  an  *  after 
the  card  number) .  Variables  which  are  vectors  are  read  sequen¬ 
tially  under  their  format  specification,  as  noted  by  the  index 
after  the  variable  name.  Special  notes  on  some  variables  are 
included  in  subsection  3.1.3  where  further  explanation  of  the 
variable  is  required.  The  format  used  in  the  following  input 
card  description  is: 

CARD  NO.,  NAME,  VECTOR  INDEX,  FORMAT,  DESCRIPTION. 
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FTIME ,  DTIME,  STIME 


(12,  8X,  2G10.5) 


FTIME  -  Final  time  of  simulation 
DTIME  -  Time  step  for  simulation 
STIME  -  Starting  time  for  simulation 
(see  note  1) 

IPRNT(I) ;  I  =  1,  9  (912) 

IPRNT  -  Print  control  vector,  print  state  variable 
IPRNT(I)  every  time  step. 

NPLTM  (12) 

NPLTM  -  Plot  control  flag,  NPLTM  =  number  of  plots 
to  be  made. 

XLAB(I);  I  =  1,  40  (40A2) 

Lable  card  for  simulation  output 
ICNTL(I);  I  =  1,  16  (1611) 

ICNTL(I)  =  Execution  option  selection 

ICNTL(l)  =  1  for  dynamic  cushion  pressure 

ICNTL ( 2 )  =  1  for  dynamic  trunk  pressure  and  fan 
characteristics 

ICNTL (3)  =  1  for  trunk  orifice  flow 

ICNTL(4)  =  1  for  separation  point  interpolation 

ICNTL (5)  =  1  for  static  trunk  shape 

ICNTL (6)  =  1  for  nonvariable  pressure  profile  (see 
card  25) 

ICNTL (1)  to  ICNTL (3)  have  to  be  1  for  complete 
simulation 


6. 


NODES,  ICFLAG,  IFXN,  IFSEP,  ILENG,  INSEP 


(1015) 


NODES  -  Number  of  trunk  nodes 

ICFLAG  -  Trunk  shape  initial  condition;  1  =  compute, 

0  =  read  initial  condition  in  terms  of  X,  Y 
(see  card  No.  12,  13) 

IFXN  -  Number  of  steps  tc  be  skipped  in  printout 

IFSEP  -  Separation  point  model  option 

=  1  for  fixed  gap  to  separation  point  height 
=  2  for  diffuser  model 

=  3  for  separation  point  fixed  at  node  number 
"INSEP" 

ILENG  -  Initial  length  of  segments  option 

=  0  compute  length  of  node-node  segments 
=  1  read  data  cards  for  the  initial  length  (see 
card  No.  11) 

INSEP  -  Separation  point  node  (IF  IFSEP  =  3) 

7.  A,  B,  L,  HYI  (8G10.5) 

A  -  Horizontal  separation  between  the  attachment  points. 
B  -  Vertical  separation  between  the  attachment  points. 

L  -  Trunk  membrane  length,  between  the  attachment  points. 
HYI  -  Initial  trunk  height. 
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8. 


SSLENG,  TPERIM 


(8G10.5) 


SSLENG  -  Length  of  trunk  with  significant  gap  area 
for  flow  to  atmosphere  from  cushion 

TPERIM  -  Perimeter  of  trunk,  around  complete  ACLS 

9.  ATC,  ATRIM  (8G1.0.5) 

ATC  -  Area  of  trunk  to  cushion  flow,  not  including  trim 
valve  (ATRIM)  (ATC  not  used  if  option  3  used) 

ATRIM  -  Area  of  trim  valve  for  fixed  flow  area  to  cushion 

10.  ATCF(I);  1=1,  nodes  +  1  (8G10.5) 

ATCF  -  Flow  area  of  trunk  orifice  per  unit  width  for 
node  to  node  link 

(assume  uniform  density  for  each  link) 

11.  *RLENGO ( I ) ;  1=1,  nodes  +  1  (8G10.5) 

RLENGO  -  Element  lengths  at  zero  extension;  needed  if 
ILENG  =  1  (see  card  No.  6) 


12.  *X ( I ) ;  1=2,  nodes  +  L  (8G10.5) 

X  Node  position  values;  needed  if  ICFLAG  =  0 
(see  card  No.  6) 

13.  *Y (I) ;  1=2,  nodes  +  1  (8G10.5) 

Y  Node  position  values;  needed  if  ICFLAG  =  0 
(see  card  No.  6) 

14.  RMASS(I);  1=1,  nodes  (8G10.5) 


RMASS  -  Nodal  mass  values 
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15.  RKVEC(I);  1=1,  nodes  (8G10.5) 

RKVEC  -  Trunk  elasticity  of  membrane  links  between  nodes 

16.  RBVEC(I);  1=1,  nodes  (8G10.5) 

RBVEC  -  Flexural  stiffness  per  node 

17.  DAMP (I);  1=1,  nodes  (8G10.5) 

DAMP  -  Nodal  damping  ratio  (see  note  2) 

18.  TREST,  DAMPR  (8G10.5) 

TREST  -  Simulation  time  damping  ratio  change 

DAMPR  -  Damping  ratio  change  multiplication  factor 
(see  note  3) 

19.  I EXT ,  RKEXTX ,  RKEXTY  (12,  2G10 . 5) 

IEXT  -  Node  of  external  spring  attachment 

RKEXTX  -  X  direction  spring  constant 

RKEXTY  -  Y  direction  spring  constant  (see  note  4) 

20.  AIFAN ,  TKVOL ,  VCH  (8G10.5) 

AIFAN  -  Inertance  of  air  in  the  fan  (see  note  5) 

TKVOL  -  Trunk  volume 

VCH  -  Cushion  volume 
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21.  CQO,  CQ1,  CQ2,  CQ3 ,  CQ4  (8G10.5) 

CQO  -  Fan  polynomial  coefficient  (see  Figure  5) 

CQ1  -  Fan  polynomial  coefficient 
CQ2  -  Fan  polynomial  coefficient 
CQ3  -  Fan  polynomial  coefficient 
CQ4  -  Fan  polynomial  coefficient  (see  note  7) 

22.  CTC ,  CTRIM,  CGAP ,  TSEP  (8G10.5) 

CTC  -  Discharge  coefficient  for  the  trunk  orifices 
CTRIM  -  Discharge  coefficient  for  the  trim  valve 
CGAP  -  Discharge  coefficient  for  the  gap 
TSEP  -  Separation  angle  (in  radians) 

23.  YGRNDS ,  SRATIO,  YDMIN  (8G10.5) 

YGRNDS  -  Hard  surface  clearance 

SRATIO  -  Constant  gap  to  separation  point  height  ratio 

YDMIN  -  Maximum  allowable  trunk  height  for  the  minimum 
gap  area  method  of  flutter  suppression. 

(strip  or  puck  induced  minimum  gap  area) 

24.  PTK ,  PCH ,  QGAP  (8G10.5) 

PTK  -  Trunk  pressure  initial  condition  (or  constant) 
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Figure  5.  Fan  pressure  versus  flow  polynomial. 


PCH  -  Cushion  pressure  initial  condition  (or  constant) 

QGAP  -  Exit  flow,  initial  condition  (not  required 
unless  ICNTL(l)  and  ICNTL(2)  are  1) 

25. *  PEXT(I);  1=1,  nodes  (8G10.5) 

PEXT  -  External  pressure  at  nodes,  used  for  initial 
shape  computation  (ICNTL(6)  has  to  be  =  1) 

26. *  NVPLOT ,  NPLT(I),  I  =  1,5  (II,  4X,  512) 

NVPLOT  -  Number  of  variables  on  plot 

NPLT  -  Variable  to  be  plotted  (state  variable  number 
to  be  plotted  =  NPLT (I)  *2+4,  see  note  6) 
(see  note  below  and  note  6) 

27. *  TPLSRT ,  TPLSTP,  DTPLOT,  XMIN,  XMAX  (6F12.5) 

TPLSRT  -  Plot  start  time 

TPLSTP  -  Plot  end  time 

DTPLOT  -  Time  increment  between  plot  points 

XMIN  -  Minimum  value  of  plot  axis  (optional) 

XMAX  -  Maximum  value  of  plot  axis  (optional) 

If  XMIN  =  XMAX  =  0.0,  plots  will  be  autoscaled  by 
program.  (See  note  below) 

*Note:  Cards  26  and  27  are  input  only  if  NPLTM  >  0 
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3.1.2  Internal  Data 


Internal  constants  used  in  the  program  include: 

CKK  -  Polytropic  expansion  exponent,  k,  1.4 

2 

G  -  Gravitational  acceleration,  g,  32.174  ft/sec 

PI  -  it,  3.1415926535 

PAT  -  Atmospheric  pressure  (absolute) ,  Pat, 

2116.8  lb/ft2 

RHO  -  Air  density,  p,  0.00234151  slugs/ft2 
3.1.3  Special  Notes  on  Input  Data 

Note  1  -  Time  step  for  simulation  has  to  be  judiciously 
chosen  by  the  user.  Using  too  large  a  time  step  will  cause 
numerical  instability,  using  too  small  a  time  step  will  make 
the  simulation  uneconomical  and  inaccurate.  In  fact,  the  time 
step  should  be  "small"  compared  to  the  smallest  period  of  system 
vibration. 

For  the  simulation  results  described  in  this  manual,  a 
time  step  of  0.001  sec  was  used.  This  may  serve  as  a  good 
starting  point  for  initial  calibration  simulations  for  a  dif¬ 
ferent  ACLS  system.  If  the  variables,  particularly  the  higher 
frequency  variables,  such  as  the  node  displacements  (X,  Y)  show 
a  rapidly  fluctuating  characteristic,  the  time  step  should  be 
reduced  until  such  tendencies  disappear.  On  the  other  hand, 
the  time  step  may  be  increased  if  the  time  step  of  0.001  sec  is 
much  smaller  than  the  time  step  at  which  the  fluctuations  appear. 
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Note  2  -  The  global  damping  values  (ratios)  specified  on 
input  must  be  chosen  by  making  an  engineering  estimate  of  the 
energy  dissipated  by  trunk  motion.  Actual  extension  and  bending 
motion  damping  values  have  been  approximated  by  a  global  (X,  Y) 
velocity  damping  force. 

Note  3  -  TREST,  DAMPR  -  TREST  is  the  simulation  time  that 
the  program  can  change  the  global  damping  ratio  values  DAMP (I) 
that  were  given  at  input.  At  time  TREST  the  values  of  DAMP (I) 
are  multiplied  by  the  constant  DAMPR.  This  option  allows  the 
initial  shape  of  the  trunk  and  fluid  flow  variables  to  be 
computed  for  a  highly  damped  system  to  get  the  required  initial 
conditions  and  then  to  simulate  the  lightly  damped  system  from 
time  equals  TREST  onward. 

Note  4  -  IEXT,  RKEXTX,  RKEXTY  -  This  option  allows  the 
attachment  of  a  spring  from  node  number  IEXT  to  the  ACLS  frame. 

The  X  and  Y  direction  stiffnesses  are  input  as  RKEXTX  and  RKEXTY. 
The  initial  trunk  shape  is  assumed  at  the  zero  extension  position. 

Note  5  -  Fan  air  inertance  is  inertance  of  the  air  residing 
in  the  fan  at  any  instant.  A  good  estimate  of  the  inertance  is 
obtained  by: 

If  =  pi  =  AIFAN 
A 


where 

1^  =  air  inertance 

=  average  flow  path  length  in  the  fan,  which  may  be 
approximated  as  the  length  of  the  fan 
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p  =  average  air  density 
A  =  Cross  section  area  of  flow. 

Note  6  -  The  state  variables  for  the  simulation  are  stored 
in  a  vector.  Each  node  requires  four  state  variables  for 
integration: 

I  =  (J  +  1)  *  4  +  1  for  node  No.  J 

STATE (I)  =  X  Velocity  of  node  (J) 

STATE (I  +1)  =  X  Displacement  of  node  (J) 

STATE (I  +2)  =  Y  Velocity  of  node  (J) 

STATE (I  +3)  =  Y  Position  of  node  (J) 

If  optional  dynamic  trunk  or  cushion  pressure  are  used  these 
variables  are  appended  at  the  node  state  vector: 

I  =  (NODES  +2)  *  4  +  ICNTL(l)  +  ICNTL(2)  *  2 
STATE (I)  =  Dynamic  cushion  pressure 

STATE (I  -  1)  =  Dynamic  trunk  pressure 

STATE (I  -  2)  =  Dynamic  fan  flow 

STATE (I  +1)  =  Trunk  membrane  length 

STATE (I  +2)  =  Flow,  trunk  to  cushion 

STATE (I  +3)  =  Flow,  trim 
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STATE (I  +  4) 


Flow,  cushion  to  atmosphere 


STATE (I  +5)  =  Average  cushion  pressure  P(t) 

STATE (I  +6)  =  Flow,  cushion  and  trunk  orifices  to 

atmosphere 

STATE (I  +7)  =  Flow,  trunk  to  atmosphere. 

Note  7  -  CQO ,  CQl,  CQ2,  CQ3 ,  CQ4  -  These  fan  flow  variables 
are  coefficients  for  a  fourth  order  polynomial.  They  are  found 
by  doing  a  linear  polynomial  regression  on  the  fan  data  using  a 
standard  regression  program.  (See  IBM-SSP  manual.) 

Note  8  -  S RATIO  -  This  ratio  is  an  optional  flow  separation 
calculation  technique  which  bases  the  height  of  the  separation 
point  gap  as  a  ratio  of  minimum  gap  to  separation  gap  equal  to 
SRATIO . 


The  last  seven  variables  are  not  state  variables  but  auxiliary 
storage  variables. 


3.1.4  Program  Option  Operation 


The  dynamic  flutter  simulation  program  includes  a  number  of 
model  options  which  allow  the  user  to  simulate  a  variety  of  ACLS 
designs  and  to  study  a  number  of  system  features.  The  simulation 
options  include: 


STATIC  PRESSURE  LOAD  PROFILE 

FLOW  INDUCED,  DYNAMIC,  PRESSURE  PROFILE 

VARIOUS  FLOW  SEPARATION  POINT  MODELS 

EXTERNAL  SPRING  ATTACHMENT 

DAMPING  RATIO  STEP  CHANGE  VERSUS  TIME 

TRUNK  SHAPE  INPUT  OR  CALCULATION 
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FLOW  DYNAMICS  INTEGRATION  WITH  STATIC  TRUNK  SHAPE 
VARIOUS  PRESSURE  DYNAMIC  MODELS 

These  features  allow  a  number  of  types  and  phases  of  analysis 
to  be  performed  with  the  program.  A  description  of  the  options 
and  how  to  select  them  is  presented  here. 

Program  options: 

1.  Static  pressure  load  profile. 

ICNTL ( 6 )  =  1 

PEXT(I)  =  Node  pressure  (external)  values.  (Optional  input) 

This  option  is  useful  for  measuring  pressurized  trunk  shapes 
without  pressure  load  dynamics  on  trunk  surface,  (for  example: 
out-of-ground-effect) . 

2.  Flow  induced  dynamic  pressure  profile. 

ICNTL (6)  =  0  (DEFAULT  MODEL) 

This  option  uses  a  Bernoulli  flow  equation  for  flow  from  the 
cushion  area  to  atmosphere.  This  is  the  normal  default  flow 
model.  When  trunk  flow  orifices  are  included  a  modification  to 
the  flow  equations  is  required: 

ICNTL (3)  =  1 

This  results  in  an  interactive  flow  computation  which  includes 
cushion  exit  flow  and  trunk  orifice  flow  which  combine  and  flow 
to  atmosphere. 

3.  Flow  separation  point  models. 

Several  models  exist  for  the  flow  separation  point  calculation. 
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a.  Diffuser  model. 

IFSEP  =  2 

This  model  uses  a  diffuser  slope  value  of  TSEP  radians  to  set 
the  separation  point  at  the  node  where  the  trunk  surface  slope 
is  closest  to  TSEP.  Typical  values  are  in  the  range  of  6  to  12 
degrees . 

b.  Fixed  gap  ratio. 

IFSEP  =  1 

This  model  sets  the  separation  point  at  the  node  point  where  the 
gap  height  is  closest  to  the  ratio  of  minimum  gap  height  divided 
by  the  ratio,  SRATIO.  (See  note  8.) 

c.  Fixed  separation  node. 

IFSEP  =  3 

INSEP  =  separation  node  number 

This  option  allows  the  user  to  set  the  separation  point  at  the 
specific  node.  This  feature  is  useful  for  simulating  strakes 
or  other  flow  separation  inducing  devices. 

d.  Trunk  orifice  flow  induced  separation. 

ICNTL { 3 )  =  1 

IFSEP  =  2 

For  the  diffuser  model  a  special  case  can  exist  when  the  trunk 
flow  orifices  are  utilized.  If  separation  would  occur  after  the 
trunk  orifice  area  it  will  actually  occur  at  the  last  trunk 
orifice  row. 
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4. 


External  spring  attachment. 


I EXT  =  node  # 

RKEXTX  =  X  direction  spring  stiffness  (linear  constant  value) 
RKEXTY  =  Y  direction  spring  stiffness  (linear  constant  value) 

This  option  is  used  to  simulate  the  attachment  of  an  external 
spring  to  the  trunk  membrane  from  the  ACLS  frame.  If  IEXT  =  0 
the  option  is  overridden. 

5.  Damping  ratio  step  change  versus  time. 

TREST  =  time  of  damping  change 

DAMPR  =  damper  ratio  multiplication  factor 

This  option  allows  the  global  damping  ratio  of  the  trunk  nodes, 

DAMP ( I ) ,  to  be  changed  during  the  simulation.  This  feature  can 
be  used  to  integrate  the  trunk  equations  in  an  overly  damped 
manner  to  reach  an  equilibrium  value  or  to  compute  quasi-static 
trunk  dynamics.  When  the  initial  conditions  of  the  system  are 
not  well  known  the  initial  computation  with  a  high  damping  ratio 
will  allow  them  to  stabilize  before  simulating  a  flutter  situation. 

6.  Trunk  shape  selection. 

ICFLAG  =  0  Input  node  x,  y  values 

=  1  Compute  node  x,  y  values 

ICFLAG  =  0  Input  node  x,  y  values 

=  1  Compute  node  x,  \  values 

The  TRUNK  subroutine  can  only  calculate  the  positions  of  equi- 
spaced  trunk  nodes.  The  user  inputted  trunk  node  x,  y  values 
are  the  most  useful  and  flexible  because  it  allows  variable 
node  placement  and  spacing.  RLENGO  is  also  controlled  in  a 
similar  manner. 
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ILENG  =  0  Compute  lengths 
=  1  Input  lengths 

The  computed  lengths  are  only  good  for  equi-spaced  nodes  and 
user  input  is  the  recommended  system. 

7.  Flow  dynamics  with  static  trunk  profile. 

ICNTL ( 5 )  =  1 

This  option  allows  flow  and  pressure  data  to  be  computed  dynami¬ 
cally,  but  for  the  trunk  shape  to  be  held  fixed.  This  feature 
is  useful  for  checking  out  flow  dynamics  with  a  controlled  static 
trunk  system. 

8.  Various  pressure  dynamic  models. 

Three  different  fluid  system  models  are  included  in  the 
program : 

ICNTL (1),  ICNTL ( 2 ) 

a.  The  pressures  can  be  fixed  to  the  input  values. 

ICNTL (1)  =  ICNTL ( 2 )  =  0 

b.  The  cushion  pressure  can  be  a  dynamic  function  of  the 
flows  and  orifice  areas  in  the  ACLS  system. 

ICNTL  (1)  =  1,  ICNTL ( 2 )  =  0 

c.  The  cushion  and  trunk  pressures  and  fan  flow  are  dyna¬ 
mic  functions  of  the  fan  dynamics,  orifice  areas,  and 
exit  flows.  If  the  dynamic  trunk  and  fan  option  is 
used,  the  dynamic  cushion  pressure  must  be  used  also. 

ICNTL (1)  =  ICNTL (2)  =  1 


28 


3.2  Program  Output 


The  printout  includes  the  following  data: 

a.  Input  Parameters  -  Trunk  parameters,  fan  parameters, 
control  parameters,  simulation  parameters,  flow  path 
parameters,  structural  parameters  are  printed  out 
after  input.  The  options  activated  by  the  ICNTL  vector 
are  printed  out  to  indicate  the  features  of  the  model 
being  simulated. 

b.  Dynamic  Simulation  Data  -  During  the  dynamic  simulation 
the  trunk  shape,  pressures  and  flow  values,  and  other 
requested  data  are  printed  out  every  IFXN  time  steps. 

c.  Time  Response  Plots  -  After  the  simulation  reaches  its 
final  time  the  program  can  produce  a  printer  plot  of  any 
of  the  dynamic  variables  such  as  trunk  node  position, 
dynamic  pressures,  or  flow  etc. 

Listing  of  output  data  in  sequential  order. 

Integration  start  time,  STIME;  integration  final  time, 

FTIME;  integration  time  step,  TSTEP. 

System  state  variables,  to  be  printed  (see  card  No.  2). 

Simulation  label,  XLAB(I). 

Control  vector,  ICNTI, (I). 

Options  in  effect  for  simulation, 

a.  Dynamic  cushion  pressure 

b.  Dynamic  trunk-fan  pressure 

c.  Trunk  orifice  flow 


d.  Separation  point  interpolation 

e.  No  trunk  motion  test 

f.  Pressure  profile. 

6.  Number  of  nodes,  NODES;  number  of  system  state  variables, 

NS TATE ;  X,  Y  node  point  selection  option  ICFLAG;  print 
out  skip  number  IFXN;  separation  point  selection  model, 

IFSEP ;  separation  point,  INSEP;  element  length  option, 

ILENG. 

7.  Horizontal  attachment  separation,  A;  vertical  attachment 
separation,  B;  trunk  length,  L;  trunk  height,  HYI . 

8.  Trunk  gap  exit  flow  length,  SSLENG;  trunk  perimeter,  TPERIM . 

9.  Trunk  to  cushion  flow  area,  ATC;  trim  valve  area  ATRIM. 

10.  Trunk  orifice  element  flow  areas,  ATCF(I). 

11.  Trunk  element  zero  extension  length,  RLENGO ( I )  (option). 

12.  X  node  position  values,  X(I). 

13.  Y  node  position  values,  Y(I). 

14.  Trunk  nodal  masses,  RMASS(I). 

15.  Trunk  elastic  stiffness,  RKVEC ( I ) . 

16.  Trunk  bending  stiffness,  RBVEC(I). 

17.  Nodal  dampina  ratios,  DAMP ( I ) . 

lb.  Damping  reset  time,  TREST,  reset  factor  DAMPR. 
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19.  External  spring  attachment  node,  IEXT;  X  direction  spring 
constant,  RKEXTX ;  Y  direction  spring  constant,  RKEXTY . 


20.  Pan  air  inertance,  AIFAN;  trunk  volume  TKVOL;  cushion 
volume  VCH . 

21.  Fan  polynomial  coefficients,  CQO,  CQ1,  CQ2,  CQ3,  CQ4 . 

22.  Trunk  to  cushion  flow  discharge  coefficient,  CTC;  trim  valve 
flow  discharge  coefficient,  CTRIM;  cushion  flow  area  dis¬ 
charge  coefficient,  CGAP ;  separation  angle,  TSEP. 

23.  Hard  surface  clearance,  YGRNDS ;  gap  to  separation  point 
height  ratio,  SRATIO;  maximum  allowable  trunk  height, 

YDMIN . 

24.  Trunk  pressure  initial  condition  or  constant,  PTK ;  cushion 
pressure  initial  condition  or  constant,  PCH ;  exit  flow  (fan) 
initial  condition,  QFAN . 

25.  External  pressure  profile,  PEXT(I)  (option). 

NOTE 

Outputs  26  to  30  are  printed  every  IFXN  time 
steps . 

26.  Simulation  time,  TIME. 

27.  Dynamic  cushion  pressure,  STATE (N ) ;  trunk  length,  STATE(N+1) 
trunk  to  cushion  flow,  QTC ;  trim  flow,  QTRIM;  cushion  to 
atmosphere  flow,  QCA.  (option,  ICNTL(l)  =  1.) 

28.  Y  node  values,  Y(I) . 
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29.  Cushion  separation  node,  ICS;  exit  separation  node,  ISEP; 
lowest  node,  INODE,  separation  point  area,  YASEP;  minimum 
gap  area,  YGAPM;  gap  area  at  separation  node,  AGAP(ISEP); 
exit  flow  velocity,  VEXIT;  exit  flow,  QEXIT;  cushion  pres¬ 
sure,  PCH;  trunk  pressure,  PTK;  fan  flow,  QFAN. 

30.  External  pressure  at  nodes,  PEXT(I). 

NOTE 

Outputs  31  and  32  are  for  plot  outputs. 

31.  Plot  curve  state  variable  numbers,  IDUM(I). 

32.  Print  plot  output  (see  sample  output). 

33.  State  variable  maximum  and  minimum  values,  YMAX(I);  YMIN(I) 
(option,  if  any  IPRNT(I)  /  0). 
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4. 


ILLUSTRATIVE  SIMULATION 


The  following  describes  the  input  data  and  the  output  print¬ 
out  and  print  plot  for  a  flutter  simulation  of  a  typical  ACLS. 

The  typical  case  simulated  includes  the  effects  of  dynamic  trunk 
and  cushion  pressure,  fan  dynamics,  and  trunk  orifice  flow.  The 
separation  point  calculation  is  based  on  the  diffuser  model  and 
no  external  spring  is  used.  The  input  variables  are  shown  in 
Figure  6  and  the  resulting  printout  is  shown  in  Figure  7. 
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Figure  6.  Input  data  test  case. 
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Figure  6.  Input  data  test  case.  (Continued) 
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3532*  1.71C  8  2. 0<«G?  2.2987  2.5"12  2.6431  2.7214  2.7324  2.6512 

•132*  1.7629  1.3357  .85894  -.*0383  -.94232  -1.2720 

Figure  7.  Program  output.  (Continued) 
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Figure  7.  Program  output.  (Continued^ 


Figure  7.  Program  output.  (Continued' 
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Figure  7.  Program  output.  (Continued) 


Figure  7.  Program  output.  (Continued) 
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Figure  7.  Program  output.  (Continues) 
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Figure  7.  Program  output.  (Continued) 


APPENDIX  A 


EQUATIONS  USED  IN  THE  MODEL 


This  program  incorporates  the  nonlinear  relationships  be¬ 
tween  the  motion  of  the  lumped  trunk  masses  and  the  system  param¬ 
eters,  such  as  the  trunk  elasticity,  damping  and  the  pressure 
forces  arising  due  to  the  fluid  flow  under  the  trunk,  taking  into 
account  the  geometry  of  the  trunk.  Figure  A-l  shows  a  schematic 
diagram  of  the  model.  As  shown  in  the  model,  the  lumped  masses 
are  connected  by  springs  representing  elasticity  of  the  trunk. 
Also  included  are  the  pressure  forces  acting  between  the  lumped 
masses  which  are  divided  between  the  adjacent  masses.  The  system 
damping  is  represented  in  this  initial  model  by  global  dampers, 
which  develop  opposing  forces  proportional  to  the  absolute  veloc¬ 
ities  of  the  lumped  masses. 

The  acceleration  components  of  a  trunk  mass  are  calculated 
along  the  X  and  Y  axes  from  summation  of  the  stiffness,  pressure 
and  damping  forces  acting  along  the  respective  axes.  Double 
integration  of  the  accelerations  gives  positions  which  are  then 
plotted  to  obtain  instantaneous  trunk  shapes. 

The  equations  used  in  the  model  are  summarized  in  the 
following : 


A . 1  Geometry  Relations  (Figure  A-2) 


h  =  Ta"'1  (<Yi+i  -  Yi)/lxi+i  -  V) 

>i  -  Tan'1  ((Vj  -  Vi+1)/(Xl  -  Xi+1))  (A-2) 


BENDING  FORCES 


Figure  A-l.  Trunk  representation  for  the  dynamic  simulation  model. 


£  . 

l 


(A- 3 ) 


i 

i 


(Yi+1  -  v2) 


Note:  Tan  1  is  a  four  quadrant  function,  therefore,  if 

6  =  Tan'1  (AY/AX),  -  it  <  6  <  tt  . 

A. 2  Force  Relations 

Spring  Force  (Figure  A-3) 


fxk.  =  coSU..)*K.  *  ,n. 
1 


+  cos (6i+i) *Ki+1  *  A£i+1 


YK. 

l 


,  ..  *  A  . 

ll  l 


+  sin<"i+1)*K.+1  *  A?'i  +  1 


where 


A  ?. . 


i 


£  . 

l 


£  . 
oi 


( A—  4 ) 


(A-  5 ) 


£  .  (£  .  is  initial  £.) 


Bending  Forces  (Figure  A-4) 


TORQUE,  T  =  K_  *  A\l> 
ti 


or 


T  =  KB*(^o  "  'P) 


(A-6 ) 


since 


T  =  F  *  x. 


where 


F  *  £  =  Kb(<Pq  - 


KB%  ~  ^ 


KB  =  K£ *  TTTT. 


so 


Ktt(»o  -  «'> 

:7i  +  H+i1 


Ki.%  -  •'■>2 


•  ,  ( p  .  +  £ .  J 

1  +  .1  1  1  +  1 


(A- 7) 


(A-8 ) 


(A-9) 
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(ALL  ANGLES  ARE  (+)  COUNTERCLOCKWISE  ROTATION) 
gure  A-4.  Bending  force  representation. 


r 


Converting  forces  and  into  component  Fx, 

at  each  node  (i)  gives  (see  Figure  A-5) 


FX.  -  -Fl  C08(6.) 


Sr 

*r 

5 


l  ' 


Fy  =  -Fl  sin^) 


F  =  -F,  cos  ( 3  .  +  tt  )  -f  cos  ( 6  •  , -i  ) 

xi+l  11  z 


Fy  =  -Fl  sin(Bi  +  tt)  -f  2  sin(Bi+1) 


FXi+2  *  ~F2  cos(‘Vl  +  "> 


FYi+2  “  'F2  sin<Si*l  +  ” 


Due  to  orthoganality : 


rr 

■'  i 


rr 

2 


3 


i+1 


IT 

i  +  1  “  2 


Fy  forces 


(A-10) 


(A-ll) 


(A-12 ) 


(A-13) 


(A-14 ) 


( A- 1 5 ) 


(A-16) 


(A-17) 
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From  trigonometric  identities: 


cos ( 9 ^  -  j)  =  sin{9i)  (A-18 ) 

sin (0^  -  j)  =  -cos  (9^)  (A-19) 

cos(6i  +  ^)  =  ~sin(0i)  ( A-20 ) 

sin(0^  +  j)  =  cos  (6i)  (A- 21 ) 


So  replacement  of  £  terms  with  equivalent  0  terms  gives: 


Fx  =  -F x  sin ( G i )  (A-22) 

Fy  =  -F±  J-cos  (0i)J  (A—  23) 


i+1 


-s in  ( 0 


i>]  - 


sin 


'  i+1 


(A-24 ) 


F 


i+1 


cos  ( 


V  ~F2 


cos (0i+1) 


( A-2  5 ) 


r 


x 


i  +  2 


-F2[-sin(ni+i) 


(A-26) 


(A-27) 


Yi  +  2 


-f2  cos(0.+1) 


Reducing  the  equations  leads  to: 


Fx  =  ~F1  sin(6i) 

Ai  " 


(A-28) 


fy  =  F1  cos(0i) 
i 


(A-29) 


Fx  =  F1  sin(0i)  -  f2  sin(0i+1) 


(A-30) 


Fy  =  -F^  costO^)  +  f2  cos(0i+1) 


(A-31) 


x  2  '  i+1' 

i  +  2 


(A-32 ) 


FYi+2  =  "F2  COS(0i+l] 


(A- 33) 


Attached  Spring  Forces 

If  an  external  spring  is  attached  to  node  I  to  suppress 
flutter,  it  creates  a  generalized  force  F’(I). 


3  1 


fg  _kext  ad 


(A- 34 ) 


In  X,  Y  components: 


FXGi  KEXT  *  AXi 


(A-35) 


F  =  -K  *  AY . 

YGi  EXT  l 


(A-36) 


where 


“AXi  =  (XGi  '  V 


(A- 37) 


-AY.  =  (Y  -  Y. ; 
1  o .  1 

i 


(A-38) 


therefore 


FXG.  KEXT  <Xo.  “  xi) 

l  l 


(A—  39) 


FYG.  =  KEXT  <Yo.  -  Yi> 


(A-40) 


Pressure  Force  (Figure  A-6) 


Defining 


P . 
l 


P  -  P 
tk  ext^ 


where  P  ^  =  Internal  Pressure  on  Membrane 

and  P  =  External  Pressure  on  Membrane  (static  or 

6X  L  • 

1 

leads  to  pressure  forces  at  node  points: 

FXPi  *  -2i[»-i*sinl6i)  +  Si  +  l*sin<0i+l>] 

fxp.  *  T  ti*=os(oi)  +  *i+1*oos(el+1)] 

1 

Damper  Forces  (Figure  A-7) 


By  definition  F  =  D  *  V 


F 


XD  . 

i 


V  *C • *  2 
x .  1 

i 


V  *£.*2 

Vi 


Vv 


K, 


5  b 


(A-41 ) 


dynamic) 


(A-42) 


( A-4  3 ) 


(A-44 ) 


(A-45) 


(A-46 ) 


Fiqure  A-7.  Damping  representation . 


.•>  s 


where 


damping  ratio  at  node  i 
nodal  mass  at  node  i 
average  stiffness  at  node  i 


£ .  = 
1 

M.  = 

l 


Differential  Equations: 


d2X . 

_ i 

dt2 


(Fxd.  +  lxk.  +  Fxp.  +  Fxb.  +  FxG.)/Mi 
11111 


d2Y . 

_ l 

dt2 


(F  ,  +  F  ,  +  F  +  F  ,  +  F  _  )/M. 

ydi  yki  yp.  ybi  i 


X . 
i 


dY  . 

_i 

dt 


Y. 

l 


A.  3  Flow  Relations 

Pressure  -  F low  Relaj  ionst.  ip  (Figure  A-8) 

Dynamic  pressure  under  membrane  modelled  by  Bernoulli 
equation  for  ideal  flow  with  no  trunk  flow. 


f>9 


(A- 4 7 ) 


(A-48) 


(A-49 ) 


flow 


P  V 

1  1 

~r  +  ~r 


P2  V2 

T  +  T" 


(A-51 ) 


where 


Q  =  V .  A .  =  V-A.,  =  V  A 

11  22  sep  sep 


(A- 52) 


therefore 


P  ,  -  P  =o 

ext^  sep  2 


f  -  Ml 


(A-53) 


at  separation  point 


P  =  P  ^  =  0 

sep  at 


P1  =  Pch'  V1  "  0 


Therefore , 


(A- 54 ) 


or 


V 

sep 


(A-55) 


s  ince 


Q 


A  *  V 
sep  sep 


(A-56) 


therefore 


(A-57) 


Pressure  Source  Dynamics 

The  system  model  includes  the  capabilities  of  fan-trunk- 
cushion  dynamics  shown  in  Figure  A-9.  The  fan  includes  a  pres¬ 
sure  versus  flow  polynomial  curve  fit  and  a  fluid  model  of  the 
trunk  and  cushion  volumes  and  orifices.  (See  final  report  for 
derivation  of  equations  A-58  through  A-61). 


D  i.  fferential  Equations  : 


d 

.it; 


( ''-FAN  ] 


P(QFAN)  ”  Ptk 


'FAN 


(A-58) 


p 


Figure  A-10.  Fan  pressure  versus  flow  polynomial. 


The  flow  model  determines  the  variations  in  pressures  and 
flows  as  a  function  of  time.  There  are  two  parts  to  the  flow 
model;  the  fluid  chambers  {that  is,  cushion  and  trunk),  and  the 
fan.  The  principal  assumptions  of  the  flow  model  are  as 
follows. 

a.  The  flow  through  all  orifices  is  one-dimensional  and 
quasi-static,  that  is,  the  pressure  in  the  plane  of 
the  orifice  is  uniform,  and  the  unsteady  state  terms 
in  Bernoulli's  equation  are  small  compared  to  the 
change  in  velocity  head. 

b.  The  flow  through  the  orifices  is  incompressible, 
that  is,  the  pressure  drop  is  small  compared  to  the 
total  pressure,  and  the  air  density  is  constant. 


M 


c.  The  pressure  and  volume  changes  of  the  air  during  ex¬ 
pansion  and  compression  in  the  various  fluid  chambers 

are  governed  by  a  polytropic  relationship,  that  is, 
k 

pv  =  const . . 

Trunk  Orifice  Flow  Effects: 

The  addition  of  flow  from  the  trunk  orifice  to  the  gap  changes 
the  pressure  profile  under  the  trunk.  The  flow/pressure  re¬ 
lations  are  iteratively  computed  between  the  cushion  and  the 
separation  point.  The  flow  is  computed  by  considering  a  number 
of  control  volumes  between  nodes.  (see  Figure  A-ll) . 

The  flow  into  each  control  volume  is  computed  as: 


qtc.  -  cTC  *  Vi .  /i[v  -  (A_62) 


'i  ',J  *  °TCi 


( A-6  3 ) 


Wi+1  =  Wi  +  dW 


(A-64 ) 


then  a  new  is  computed 


ri+l 


—  *  v. 

2  l 


X 

/  ¥ 


(A-65 ) 


The  flow-pressure  pattern  is  computed  under  the  trunk  step-by- 
step  out  to  the  separation  point.  The  flow  computation  is 
iterated  by  varying  W ( 1 )  until  the  exit  pressure  at  the  sepa¬ 
ration  point  is  within  bounds. 
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APPENDIX  B 


THE  EIGENVALUE  PROGRAM 


This  program  was  developed  as  a  part  of  the  contract.  By 
generating  eigenvalues  and  eigenvectors  of  the  trunk  for  various 
models,  this  program  can  assist  in  understanding  the  trunk 
behavior . 


B.l  Program  Description 

The  eigenvalue  program  has  three  parts.  The  first  part 
analyzes  the  two  dimensional  vibration  of  a  membrane  using 
coupled  longitudinal  and  lateral  motions  including  tension  and 
elasticity  effects.  The  second  part  analyzes  vibration  of  a 
membrane  which  does  not  have  longitudinal  motion  and  vibrates 
only  in  the  lateral  direction  similar  to  a  stretched  string. 

The  third  part  of  the  program  is  relevant  for  an  elastic  membrane 
which  can  vibrate  only  in  the  longitudinal  direction.  This  mode 
of  vibration  is  similar  to  that  of  a  bar.  The  program  can  work 
with  damped  or  undamped  systems.  In  addition,  the  program  calcu¬ 
lates  the  natural  frequencies  for  an  equivalent  string.  The  out¬ 
put  of  the  program  consists  of  the  eigenvalues  and  eigenvectors 
(optional) . 

A  list  of  the  subroutines  used  is  presented  in  Table  B-l. 

The  input  data  description  for  the  program  and  a  sample  output 
are  also  presented  in  this  Appendix.  However,  first  the  various 
models  used  in  the  program  are  described  in  Table  B-l. 


TABLE  B-l 


A  SUMMARY  OF  EIGENVALUE  PROGRAM  SUBROUTINES 


No . 

Subroutine 

Pi  lmary  Function 

Group 

1 

FMAEVEC 

Main  program;  I/O  control, 
coordinate  analysis;  form 
matrices 

*MAIN* 

2 

CLEAR 

Clear  matrix  to  zero 

MATRIX 

3 

r”TMAT 

Print  matrix 

I/O 

4 

TRUNK 

Computer  trunk  shape 

GEOMETRY 

5 

ELEMK 

Form  element  stiffness  matrix 

GEOMETRY 

6 

MOVE 

Copy  matrix 

.  . 

MATRIX 

7 

EIGPAC 

Eigenvalue  'eigenvector  computa¬ 
tion  coordination  module 

EIGEN 

8 

CMINV 

Complex  matrix  inversion 

MATRIX 

9 

PUTEIG 

Print  eigenvalues 

I/O 

10 

HSBn 

SbF  eigenvalue  routine 

EIGEN 

LI 

ATEIG 

SSP  eigenvalue  routine 

EIGEN 

12 

VECPAC 

Eigenvector  computation  and  out¬ 
put  coordination  module 

EIGEN 

13 

EVECTR 

Solve  complex  system  of  equations 

EIGEN 

B.2  Models  Used  in  the  Eigenvalue  Program 


B.2.1  Lateral  Vibration  Model 


In  Figure  B-l,  the  force  due  to  displacement 


Y 

P 


is : 


F  =  -  T  Sin  (a,)  +  T  Sin  (a  ) 

it  H  r* 


a  first  order  approximation, 


Sin  (a  ^ 


Sin  (cxp) 
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E 


(B-l) 


HARD 
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Figure  B-l. 


Lateral  vibration  model. 
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for  l 
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Then , 


(yP+l  +  yp-l} 


0 


the  boundary  conditions  are: 


y 


o 


yn  +  l  = 


0. 


Writing  equation  (B-4)  in  a  Matrix  form: 


Y  +  AY  =  0 


h 

2  A  i ,  A  j  / 
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y3 
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_  yn  _ 

Solution  det  | i |  =  0  gives  eigenvalues  of  the  vibration  modes. 


A  . 
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2 
w  . 
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The  addition  of  dampinq  to  the  model  requires  the  addition 
of  a  damper  force  it  each  node. 


where 


FDp  =  -Vp*Bp  (see  Figure  B-2) 


Bp  =  Damping  constant  of  point  P 
Vp  =  Velocity  of  point  P 


so  equation  (B-4)  becomes: 
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Figure  B-2.  Nodal  dampers  (lateral) 


assumption) 


Solution  of  det  j  >,  j  =  0  gives  the  eigenvalues  (natural 

frequencies)  of  the  vibration  modes. 


B.2.2  Longitudinal  Vibration  Model 


In  Figure  B-3  the  force  due  to  displacement  X^  is: 


F  . 
1 


-K.X. 
1  1 


K  .  ,  X  . 
l  +  l  i 


(Ki  +  Ki+1! 


Ki+lXi+l 


Ki+lXi+l 


KiXi-l  (B“9) 


K.X.  , 
l  1-1 


(B-10) 


F  =  MA 


Figure  B- 1 . 


Longitudinal  vihration  model. 


Solution  of  det  |_Xj  =  0  gives  the  eigenvalues  of  the 

vibration  modes. 


\  .  =  10  . 

1  l 

Addition  of  damping  to  the  model  requires  the  addition  of  a 
damper  force  at  each  node. 

F  .  =  -V . *B .  ( B— 1 4 

Di  li 

(see  Figure  B-4 ) . 
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Equation 


Ki 


X.,  X. 

i  1 


A/W 


m . 

j 

f 

l 

B, 


Figure  B-4.  Nodal  dampers  (longitudinal). 


( B-12 )  becomes 


B  . 

+  —  X  .  =  0 

m .  i 
l 
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In  matrix  form: 


B.2.3  Coupled  Longitudinal  and  Lateral  Motion 

The  two  dimensional  matrix  formulation  requires  a  two  co¬ 
ordinate  vector  (X  ,  Y,,)  at  ”ach  node  creating  a  2*N0DES  state 
space.  The  matrix  formulation  of  the  model  requires  lineariza¬ 
tion  of  the  equations  about  some  configuration  of  the  membrane. 
Each  spring  element  has  its  linear  stiffness  matrix  converted  to 
the  global  coordinate  frame  as  shown  in  Figure  B-5. 


Stiffness  of  an  element 
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in  one  dimension 
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in  two  dimensions 


R 


[K]i(6}i 


UT  =  [ujVj  U2v2] 

N+l 

K  =  K±  ;  KU  =  R  (B-17) 

1 


K  =  global  stiffness  matrix 
U  =  displacement  vector 
R  =  forcing  load  vector. 

The  membrane  tension  stiffness  effects  must  be  added  to  the 
element  spring  stiffness  matrix.  The  computation  of  the  two  di¬ 
mensional  equivalent  stiffness  for  tensile  forces  requires  a 
linearization  of  transverse  nodal  motion  into  the  second  coordi¬ 
nate  frame. 

From  lateral  one  dimensional  development  earlier,  effective 
lateral  stiffness  due  to  tension: 

k  =  | - y  in  Figure  B-6  (B-18) 

1 '  i  '  i-l J 


since 


k  *  cos  ( ■* )  ^  k 
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."■] 


ft  =  K(yi+i  "  2yi  +  yi-i} 

* 

t, 

* 

(y  is  rotated  relative  to  y  by  a.) 

The  lateral  y  displacement  must  be  transformed  into  (x,Y) 
frame  as  follows 

Transverse  displacement  D  transformation  into  (X,Y)  frame 

t, y  cos(a^)  +  Ax  sin(n^)  =  D  (Figure  B-7)  (B-20) 


Fiuuro  13-7. 


Force  transformation. 


y±  =  yi  Jcos  (ai)j  +  xi  [sin  (a±)j 


X . 
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=  sin  a. 

2  a 


y^ 

[sin  ai_1  x._1  +  cos  ai_1  y..^ 


(B-21) 
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(B-22) 


B.3  Eigenvalue  Input  Data 


1.  ICNTL(I),  I  =  1,  10  IPRNT(I),  1=1,  10  (2011) 


Program  Control  vectors 


ICNTL  (N) 

= 

1, 

Eigenvalues ; 

N 

= 

(1) 

compute 

N 

= 

(2) 

compute 

N 

= 

(3) 

compute 

N 

= 

(4) 

compute 

N 

= 

(5) 

compute 

N 

= 

(6) 

compute 

IPRNT  N 

1/ 

print  matrix 

Values 

2,  Eigenvectors;  0,  skip 
two-dimensional  Eigenvalues 
two-dimensional  damped  Eigenvalues 
transverse  string  undamped  Eigenvalues 
transverse  damped  Eigenvalues 
longitudinal  bar  undamped  Eigenvalues 
longitudinal  bar  damped  Eigenvalues 

as  above;  0,  No  print 


S4 


2. 


NODES 


(12) 


NODES  -  Number  of  nodes 

3.  LS,  RL,  AX,  BX,  TENSN ,  HY  (8G10.5) 

LS  -  membrane  length,  stretched 

RL  -  membrane  length,  unstretched 

AX  -  horizontal  distance  between  attachment  point 
BX  -  vertical  distance  between  attachment  point 
TENSN  -  membrane  preload  tension 

HY  -  trunk  height 

4.  MASS (I),  1=1,  nodes  (8G10.5) 

MASS  -  lumped  parameter  nodal  mass 

5.  RKVEC(I),  1=1,  nodes  +  1  (8G10.5) 

RKVEC  -  elastic  stiffness  of  element 

6.  DAMP ( I ) ,  I  =  1,  nodes  (8G10.5) 

D7\MP  -  nodal  damping  in  both  x  and  y  directions 


7 .  IXF 


(2011) 


IXF  -  element  length  select  flag 
(0,  default;  1,  read  card  8) 

(Option)  RLLNGO'i) ,  I  -  1,  nodes  +  1  (8G10.5) 

RLENGO  -  element  unstretched  length 


85 


9. 


IXY 


(2011) 


IXY  -  coordinate  point  select  flag 
(0,  read  cards  10,  11;  1,  compute) 

10.  (Option)  X(I),  1=1,  nodes  (8G10.5) 

X  -  X  values  of  mass  nodes 

11.  (Option)  Y(I),  1=1,  nodes  (8G10.5) 

Y  -  Y  values  of  mass  nodes. 

B . 4  Eigenvalue  Program  Output 

The  printout  includes  the  following  data: 

a.  Input  Parameters  -  Trunk  parameters,  structural  param¬ 

eters,  and  program  control  data  are  printed  out  after 
input.  A  sample  input  is  shown  in  Figure  B-8. 

b.  Matrix  Model  Data  (Optional)  -  The  matrices  generated 
for  the  models  can  be  printed  out.  These  matrices  are 
the  global  stiffness  matrices  of  the  algebraic  models 
for  the  trunk. 

c.  Eigenvalues  -  The  Eigenvalues  and  natural  frequencies 
of  the  matrix  models  are  printed  out  in  radian  and 
Hertz  frequencies,  respectively. 

d.  Eigenvectors  (Optional)  -  The  Eigenvectors  for  each 
Eigenvalue  of  complex  pair  of  Eigenvalues  are  printed 
and  then  the  normalized  displacement  Eigenvectors 
(velocity  terms  for  damped  models  are  neglected)  are 
printed . 


Ficjure  B-8.  Sample  input  data. 


A  list  of  program  output  variables  in  sequential  order 

1.  Number  of  nodes,  NODES 

2.  Membrane  length,  LS;  rest  length,  RL,  X  end  point,  AX; 
end  point,  BX;  membrane  tension,  TENSN 

3.  Mass  Modes,  MASS  (I) 

4.  Elastic  stiffness,  RKVEC(I) 

5.  Damping  ratio  (global)  at  node,  DAMP (I) 

t .  Membrane  element  unstretched  length,  RLENGO 

7.  Node  coordinate  positions,  X(I),  Y(I) 

8.  Continuous  string  frequencies,  W(I) 

9.  2D  stiffness  matrix,  G(I,J)  (Option) 

10.  Eigenvectors  of  matrix,  Z,  W(I),  Y(I),  Hertz,  radians/ 
sec,  eigenvalue  (real,  imaginary) 

LI.  Eigenvectors,  X(I),  Y(I)  eigenvector  (real,  imaginary) 

(Option) 

12.  Eigenvectors,  X(I),  normalized  displacement  terms  only 

(Option ) 

21  >  damped  st.;  '  incus  matrix,  A  (I ,  J)  (Option) 

4.  As  10  above,  eigenvalue 

l.'- .  A.  1 1.  above,  eigenvector 


16. 

As 

12 

above , 

normalized  eigenvector 

17. 

ID 

lateral  stiffness  matrix,  GS ( I,J)  (Option) 

18. 

As 

10 

above , 

eigenvalue 

19. 

As 

11 

above. 

eigenvector 

20. 

As 

12 

above , 

normalized  eigenvector 

21. 

ID 

damped  lateral  stiffness  matrix,  XW(I,J)  (Option) 

22. 

As 

10 

above. 

eigenvalues 

23. 

As 

11 

above , 

eigenvectors 

24. 

As 

12 

above , 

normalized  eigenvectors 

25. 

ID 

longitudinal  stiffness  matrix,  GS(I,J)  (Option) 

26. 

As 

10 

above , 

eigenvalues 

27. 

As 

11 

above , 

eigenvectors 

28. 

As 

12 

above , 

normalized  eigenvectors 

29. 

ID 

damped  longitudinal  stiffness  matrix,  XW(I,J)  (Option) 

30. 

As 

10 

above , 

eigenvalues 

31 . 

As 

11 

above. 

eigenvectors 

32. 

As 

12 

above. 

normalized  eigenvectors. 

Figure  B-9  shows  a  sample  program  output. 
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Figure  B-9.  Sample  output  data.  (Continued) 


r 


i 


B . 5  Principal  Nomenclature 


Program  Variable 
Name 

Symbol 

Explanation 

A,  XW,  XZ(I,J) 

Matrix,  global  stiffness 

AS,  G,  GS (I , J) 

Matrix,  global  stiffness 

AT,  WORKM,  XT ( I , J) 

Matrix,  global  stiffness 

AX 

a 

Horizontal  distance  between 
attachment  points 

BX 

b 

Vertical  distance  between  attach¬ 
ment  points 

COSALP (I) 

cos (8.^) 

Cosine  of  ALPHA (I) 

COSTHE ( I ) 

cos (9^) 

Cosine  of  THETA (I) 

DAMP ( I ) 

Di 

Damping  ratio  of  node  (I) 

ESTIF (I, J) 

Matrix,  element  stiffness 

EVEC (I) 

Work  vector  space 

ICNTL(I) 

Program  control  vector 

IPRNT(I) 

Printer  control  vector 

IXF 

Flag,  element  initial  length 
selection 

IXY 

Flag,  X,  Y  position  selection 

LS 

l 

s 

Trunk  length 

LI 

h 

Trunk  inner  length 

L2 

l2 

Trunk  outer  length 

MASS (I) 

mi 

Nodal  mass 

NODES 

Number  of  nodes 

PHI1 

♦l 

Trunk  angle,  inner 

PHI2 

02 

Trunk  angle,  outer 
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Program  Variable 
Name 

Symbol 

R1 

R1 

R2 

R2 

RKVEC(I) 

RK. 

1 

RL 

RLENG(I) 

RL. 

RLENGO ( I ) 

RL0i 

SINALP (I) 

sinfe^ 

SINTHE(I) 

sin  (8^) 

TENSN 

T 

W 

w 

W0RK1 (I) 

W0RK2 (I) 

X(I) 

Xi 

Y  ( I ) 

Y. 

1 

Explanation 
Trunk  radius,  inner 
Trunk  radius,  outer 
Trunk  elasticity,  element  (I) 

Trunk  element  length  (I) 

Trunk  element  initial  length  (I) 
Sine  of  ALPHA (I) 

Sine  of  THETA (I) 

Tension  in  trunk 
String  frequency 
Matrix,  work  space 
Matrix,  work  space 
X  position  of  node  (I) 

Y  position  of  node  (I) 
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APPENDIX  C 


PRINCIPLE  PROGRAM  NOMENCLATURE 

The  variables  used  in  the  flutter  simulation  program  are 
defined  in  this  appendix.  Also  mentioned,  corresponding  to  the 
appropriate  computer  program  variables,  are  the  symbols  used  in 
the  analysis  of  the  trunk  model.  All  program  variables  are  in 
ft- lb-sec  units  except  where  indicated  to  the  contrary. 


Program  Variable 


Name 

Symbol 

Explanation 

A 

a 

Horizontal  distance  between  inner 
and  outer  trunk  attachment  point 

AGAP (I) 

Ai 

Flow  gap  trunk  to  ground 

A  IF  AN 

Xf 

Fan  inertance 

ATC 

Atc 

Area  trunk  to  cushion 

ATCF (I) 

htcfi 

Area  trunk  to  cushion  element 

ATRIM 

Atr 

Area  trim  valve 

B 

b 

Vertical  distance  between  inner 
and  outer  trunk  attachment  point 

CGAP 

cg 

Discharge  coefficient  gap  flow 

CKK 

Polytropic  expansion  coefficient 

COSPHI (I) 

Cos  {(})i) 

Cosine  of  PHI  (I) 

COSTHE (I) 

Cos ( 0 i ) 

Sine  of  THETA  (I) 

CQO 

ao 

Fan  polynomial  coefficient  1 

CQl 

al 

Fan  polynomial  coefficient  2 

CQ2 

Olj 

Fan  polynomial  coefficient  3 

CQ3 

u3 

Fan  polynomial  coefficient  4 

CQ4 

a4 

Fan  polynomial  coefficient  5 

Program  Variable 

Name  Symbol 


Explanation 


CTC 

Ctc 

Discharge  coefficient  trunk  to 
cushion,  flow  other  than  trim. 

CTRIM 

Ctr 

Discharge  coefficient  trim  value 

DAMP (I) 

Di 

Damping  ratio  X,  Y  at  node  (I) 

DAMPR 

Damper  rest  value 

DERY (I) 

it  <si> 

Derivatives  of  state  variables 

DQ 

dq 

Incremental  flow  trunk-channel 

DT1ME 

dt 

Integration  time  step 

DVCH 

St  (VCh> 

Cushion  volume  rate 

of  change 

DW  ( I ) 

dw : 

Incremental  mass  flow  trunk  - 
channel  element  (I) 

FORCXB(I) 

Fxbi 

Force,  X  direction. 

bending  (I) 

FORCXD ( I ) 

Fxdi 

Force,  X  direction. 

damping  (I) 

FORCXK (I) 

Fxk . 

1 

Force,  X  direction, 

elasticity  (I) 

FORCXP (I) 

FxPi 

Force,  X  direction, 

pressure  (I) 

FORCYB ( I ) 

Fybi 

Force,  Y  direction. 

bending  (I) 

FORCYD ( I ) 

Fydi 

Force,  Y  direction. 

damping  (I) 

FORCYK(I) 

Fyki 

Force,  Y  direction. 

elasticity  (I) 

FORCYP ( I ) 

Fypi 

Force,  Y  direction, 

pressure  (I) 

GAPMIN 

Minimum  gap  area 

HY 

Trunk  height 

ICFLAG 

Flag,  trunk  shape  selection 
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Program  Variable 


Name 

Symbol 

Explanation 

ICNTL (I) 

Control  vector  for  options 

ICS 

Cushion  separation  node 

IFSEP 

Flag,  separation  point  selection 

ILENG 

Flag,  segment  length  selection 

INODE 

Node  number  of  lowest  trunk  point 

INSEP 

Stake  location  node 

ISEP 

Separation  point  node 

L 

a 

Trunk  membrane  length 

LI 

*1 

Inner  trunk  length  (inner  attach¬ 
ment  point  to  bottom) 

L2 

Outer  trunk  length  (outer  attach¬ 
ment  point  to  bottom)  (£. ,  £_  used 
only  by  subroutine  TRUNK! 

NODES 

Number  of  mass  nodes 

PAT 

Pa 

Atmospheric  pressure 

PCH 

Pch 

Cushion  pressure 

PCRIT 

Critical  flow  pressure 

PEXT (I) 

External  pressure  on  trunk 

PHI  (I) 

♦i 

Angle  PHI  at  node  (I) 

PLOST 

Pressure  loss  due  to  momentum 

PRESUR ( I ) 

Pressure  differential  on  trunk 
element 

PTK 

Ptk 

Trunk  pressure 

QCA 

qca 

Flow,  cushion  to  atmosphere 

QEXIT 

qexit 

Flow  at  exit  (separation) 

QFAN 

qfan 

Fan  flow 

QGAP 

qgap 

Gap  flow 

QIN 

qin 

Flow  to  trunk 
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Program  Variable 
Name 


! 


Symbol  Explanation 


QINT 

Momentum  change  integral 

QOUT 

^out 

Flow,  out  of  trunk 

QTA 

qta 

Flow,  trunk  to  atmosphere 

QTC 

qtc 

Flow,  trunk  to  cushion 

QTOT 

Total  gap  flow 

QTRIM 

qtr 

Flow,  trim  valve 

RBVEC ( I ) 

RBj. 

Bending  stiffness,  node  (I) 

RHO 

P 

Air  density 

RITER 

Iteration  parameter 

RKEXTX 

External  spring  stiffness,  X 

RKEXTY 

External  spring  stiffness,  Y 

RKVEC ( I ) 

RK. 

1 

Trunk  elasticity 

RLENG  (I) 

RLi 

Trunk  element  length 

RLENGO ( I ) 

RLo. 

Trunk  element  unstretched  length 

RMASS  (I) 

X 

M. 

1 

Nodal  mass 

SIE  (I) 

V  . 

1 

Angle  SIE 

SINPHI  (I) 

SIN  (4>i) 

Sine  of  PHI (I) 

SINTHE  (I) 

SIN  (9i) 

Sine  of  THETA  (I) 

SSLENG 

Trunk  gap  length 

STATE ( I ) 

Si 

State  variable  (I) 

TEMPAT 

Atmospheric  temperature 

THETA (I) 

0 

1 

Angle  THETA (I) 

TIME 

t 

Simulation  time 

TKVOL 

vtk 

Trunk  volume 

Damper  reset  time 


TREST 


Program  Variable 


Name 

Symbol 

Explanation 

TSEP 

Separation  point  angle 

VCH 

< 

o 

p- 

Volume  of  cushion 

VEL(I) 

V. 

1 

Velocity  of  flow  at  node  (I) 

VEXIT 

V 

exit 

Initial  estimate  of  flow  velocity 
at  separation  point 

VSEP 

V 

s 

Velocity  at  separation  point 

W(I) 

w. 

1 

Mass  flow  at  node  (I) 

X(I) 

X  . 

1 

X  position  of  node  (I) 

XEXTO 

X  position  of  spring  at  rest 

XZETA 

^x 

X  damping  ratio 

Y  ( I ) 

Y. 

l 

Y  position  of  node  (I) 

YASEP 

Ys 

Flow  separation  point  area 

YCSEP 

Cushion  separation  point  area 

YDMIN 

Minimum  trunk  Y  displacement 

YEXTO 

Y  position  of  spring  at  rest 

YGAPM 

Minimum  gap  area 

YGRNDS 

Hard  surface  clearance 

YSEPX 

Separation  point  gap 
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APPENDIX  D 


i 


PROGRAM  LISTINGS 


All  programs  and  subroutines  in  this  report  have  been  de¬ 
signed  to  work  under  ANSI. 66  FORTRAN  IV  and  supersets  of  the 
former . 

•  Dynamic  Program  Listings  (subsection  D.l) 

•  Eigenvalue  Program  Listings  (subsection  D.2) 


NOTE 

Subroutines  HSBG,  and  ATEIG  have  been 
omitted.  Information  on  these  routines 
can  be  found  in  the  IBM  SSP  manual. 


| 


j 

i 

t 


D.l  Dynamic  Simulation  Programs 

The  following  programs  and  subroutines  are  included. 

Programs  -  DYSYS 
Subroutines  -  RKDIF 
EQSIM 
TRUNK 
PUTVEC 
ERROR 
PLOTTER 
PACKER 
PRNTPLOT 
PLOT 
P STORE 
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D. 2  Eigenvalue  Analysis  Program 

The  following  programs  and  subroutines  are  included. 
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